//----------------------------------------------------------------------计算几何

#define MAXP 100
const double eps=1e-8,inf=1e18,pi=3.14159265358979324;

struct point{double x,y;};
struct line{point p,q;double a;};
struct polygon{int n;point p[MAXP];};
inline int dcmp(double x){return x<eps?-1:(x>eps?1:0);}

//叉积:OA X OB
double cross_product(point &o,point &a,point &b)
{
	return (b.y-o.y)*(a.x-o.x)-(a.y-o.y)*(b.x-o.x);
}
//点积:OA * OB
double dot_product(point &o,point &a,point &b)
{
	return (a.x-o.x)*(b.x-o.x)+(a.y-o.y)*(b.y-o.y);
}
//点到直线距离:返回O到直线AB的距离,h为垂足
double pointtoline(point &o,point &a,point &b,point &h)
{
	double d=dis(a,b),
	       s=cross_product(a,b,o)/d;
	h.x=o.x+s*(b.y-a.y)/d;
	h.y=o.y-s*(b.x-a.x)/d;
	return fabs(s);
}
//两点间距离
double dis(point &a,point &b)
{
	double x=a.x-b.x,y=a.y-b.y;
	return sqrt(x*x+y*y);
}
//两点距离平方
double dis2(point &a,point &b)
{
	double x=a.x-b.x,y=a.y-b.y;
	return x*x+y*y;
}
//直线与直线的交点
point line_intersection(line &l1,line &l2)
{
	double s1=cross_product(l1.p,l1.q,l2.q),
	       s2=cross_product(l1.p,l2.p,l1.q);
	point ret;
	ret.x=(s2*l2.q.x+s1*l2.p.x)/(s1+s2);
	ret.y=(s2*l2.q.y+s1*l2.p.y)/(s1+s2);
	return ret;
}
//线段与圆的交点
int segment_circle(point &a,point &b,point &o,double r,point ret[2])
{
	if(dis2(a,o)<=r*r&&dis2(b,o)<=r*r)return 0;
	point chuiz;double tmp,h,m;int n=0;
	h=pointtoline(o,a,b,chuiz);
	if(h>=r)return 0;
	tmp=dis(a,b);m=sqrt(r*r-h*h);
	ret[0].x=chuiz.x+(a.x-b.x)*m/tmp;
	ret[0].y=chuiz.y+(a.y-b.y)*m/tmp;
	ret[1].x=2*chuiz.x-ret[0].x;
	ret[1].y=2*chuiz.y-ret[0].y;
	if(dot_product(ret[0],a,b)<0)n++;
	else ret[0]=ret[1];
	if(dot_product(ret[n],a,b)<0)n++;
	return n;
}
//三角形OAB和圆O的交的面积
double triangle_circle(point &o,point &a,point &b,double r)
{
	point p[4];int n;double ret=0,d[4];
	p[0]=a;
	n=segment_circle(a,b,o,r,p+1);
	p[n+1]=b;
	for(int i=0;i<=n+1;i++)d[i]=dis(p[i],o)-r;
	for(int i=0;i<=n;i++)
	{
		if(dcmp(d[i])>0||dcmp(d[i+1])>0)
			ret+=r*r*asin(cross_product(o,p[i],p[i+1])/((d[i]+r)*(d[i+1]+r)));
		else ret+=cross_product(o,p[i],p[i+1]);
	}
	return ret*0.5;
}

/******************半平面交:开始**********************/
bool cmp_a(line l1,line l2)//按极角排序
{
	int ret=dcmp(l1.a-l2.a);
	if(ret<0)return true;
	if(ret>0)return false;
	if(cross_product(l1.p,l1.q,l2.p)>0)return false;
	return true;
}

bool point_in_halfplane(point p,line &l)
{
	if(cross_product(l.p,l.q,p)>0)return true;
	return false;
}

void standardlize(line l[],int &n)
{
	int i,j;
	sort(l,l+n,cmp_a);
	for(i=0;i<n;i++)if(i && !dcmp(l[i].a-l[i-1].a))l[i].a=10086;
	for(i=1;i<n&&l[i].a!=10086;i++);
	for(j=i;j<n;j++)
		if(l[j].a!=10086)l[i++]=l[j];
	n=i;
}

void half_plane_intersection(line l[],int &n,polygon &ret)
{
	int i,j,k,top,bot;
	standardlize(l,n);//极角排序,去重
	top=2;bot=1;
	dq[top]=l[1];
	dq[bot]=l[0];
	for(i=2;i<n;i++)
	{
		while( 	top>bot &&
		       	!point_in_halfplane(line_intersection(dq[top],dq[top-1]),l[i]))
			top--;
		while( 	top>bot &&
		       	!point_in_halfplane(line_intersection(dq[bot],dq[bot+1]),l[i]))
			bot++;
		dq[++top]=l[i];
	}
	j=top+1;k=bot+1;
	while(j!=top||k!=bot)
	{
		j=top;k=bot;
		while(	top>bot &&
			!point_in_halfplane(line_intersection(dq[top],dq[top-1]),dq[bot]))
			top--;
		while(	top>bot &&
			!point_in_halfplane(line_intersection(dq[bot],dq[bot+1]),dq[top]))
			bot++;
	}
	dq[--bot]=dq[top];
	for(i=bot,j=0;i<top;i++,j++)ret.p[j]=line_intersection(dq[i],dq[i+1]);
	ret.n=j;
}
/***********************半平面交:结束**********************************/
